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Abstract. The plane-wave decomposition method (PWDM), a widely used 
means of numerically finding eigenstates of the Helmholtz equation in billiard 
systems is described as a variant of the mathematically well-established boundary 
integral method (BIM). A new unified framework encompassing the two methods 
is discussed. Furthermore, a third numerical method, which we call the Gauge 
Freedom Method (GFM) is derived from the BIM equations. This opens the way 
to further improvements in eigenstate search techniques. 



1. Introduction 

Solving the Helmholtz equation within a domain given Dirichlct boundary conditions 
is of great interest to both physicists 1^ and engineers. Firstly, the Helmholtz equation 
is the simplest example of a wave equation. Furthermore, this equation may be used 
to describe acoustics waves, microwave systems, and in particular the wavefunction 
of a quantal particle inside nano-scale devices [2] such as quantum-dots, where the 
motion of the electrons can be regarded as a free motion within a box. For this reason 
it has become a prototype problem in studies of quantum chaos. 

Of particular interest are the wavefunctions ^{x) of a stationary particle in a 
two-dimensional box (a so-called billiard system). These wavefunctions are solutions 
of the homogeneous Helmholtz equation H'if{x) — 0, where the differential operator 
H is defined as 

n = -V^ - fc2. (1) 

Note that for the special case k — 0, the Helmholtz equation reduces to Laplace's 
equation. Given a closed boundary we can ask whether this equation has a non-trivial 
solution that satisfies Dirichlet boundary conditions ^{x) = 0. 

Two main numerical strategies have been suggested to date in the literature 
in order to find the eigenstates of the Helmholtz equation (for more comprehensive 
reviews and references, see for example The first strategy can be described 

as a 'Laplacian diagonalization'. A basis is selected such that the functions it 
contains satisfy the Dirichlet boundary conditions. For example, in some cases 
one can use conformal mapping to determine a basis |B] (and see also [7]). The 
Laplacian operator is then written in this basis and diagonalized. Numerically, some 
truncation is required, and the diagonalization only determines all the eigenstates up 
to some maximum wavenumber k^^^. Thus, the Laplacian diagonalization strategy 
is inherently limited, and can not be used for the purpose of finding high-lying 
eigenstates. 
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The second numerical strategy, which is the object of this paper, can be described 
as a 'boundary approach'. This strategy is based on the observation that the 
eigenfunctions are completely determined by their behavior at the boundary. The 
boundary methods use basis functions that satisfy the Helniholtz equation inside the 
billiard at fixed k. A linear combination of the basis functions is then selected such that 
the boundary conditions are satisfied. Thus, in order to find the eigenstates, one only 
needs to study the small k window that contains the energy range of interest. Therefore 
the method is naturally suitable for the purpose of finding high-lying eigenstates. 
For 2D billiards, the Laplacian diagonalization requires 2D grid calculations. This 
is a heavy numerical task. The boundary approach, on the other hand reduces the 
calculations to a ID boundary grid. 

In the quantum chaos community, two boundary methods are commonly 
employed. The first one is referred to as the boundary integral method (BIM) |H], 
while the other is what we call here the decomposition method (DEM) , of which the 
plane-wave decomposition method (PWDM) |^ is a special case. Extensions of the 
standard PWDM have been used in and in [5] . 

Usually, the BIM and the PWDM are considered to be two independent self- 
contained procedures. Several studies have been done in order to compare their 
capabilities ^21 • While the BIM equation is exact, its convergence is very slow (power 
law in the number b of discretization points per half- wavelength) . On the other hand 
while the PWDM is mathematically limited (e.g. the maximal b is semiclassically 
determined), it is still found to be extremely efficient in practice. Hence there is 
definitely a need to develop hybrid boundary methods. 

In the present paper, we adopt a new point of view through which we regard the 
BIM and the DEM as sequences of four independent steps. By doing so, we are going 
to make the observation that the DEM and the BIM are strongly related: The two 
procedures are based on the diagonalization of literally the same matrix! As a bridge 
between them, we will highlight an intermediate strategy which we call the gauge- 
freedom method (GFM). In a follow-up paper, this framework will lead the way to 
improved eigenstate search techniques combining the strengths of the two boundary 
methods [T^ . 

Our unified description of the different boundary methods can be summarized by 
the following set of four steps that are common to the BIM and the DEM, and as we 
will show later, to the GFM: 

• Choice of a set of basis functions Fj{x; k). 

• Definition of the Fredholm matrix Ajs{k). 

• Procedure for construction of the wavefunction ^P^. 

• Definition of the quantization measure S{k). 

The first step consists of selecting a set of basis functions Fj {x; k) labeled j — 1..N. 
All boundary methods rely on basis functions that satisfy the Helmholtz equation 
inside the billiard. Thus, a superposition of such basis functions is an eigenfunction if 
it vanishes along the boundary. The choices of bases that correspond to the PWDM, 
to the primitive version of the BIM and to the simplest variation of the GFM are as 
follows: 



F,{x;k) 
F,{x;k) 
F,{x;k) 



cos(0j -|- knj ■ x) 
Yo{k\x - Xj\) 
Jo{k\x - Xj\) 



YO-BIM 



JO-GFM 



PWDM 



(2) 
(3) 
(4) 
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For the purpose of the numerical treatment we represent the boundary by a set of 
points Xs with s = 1 • ■ • M. In practice, we choose a set of M equaUy spaced points, so 
that the spacing is As = L/AI where L is the perimeter of the biUiard. Depending on 
details of the numerical strategy, the number of points along the boundary is either 
taken to be equal to the number of basis functions (Af — N), or it may be larger 
(M > N). The Fredholm matrix is defined as 

A,,(/c) = F^{x,;k) (5) 

Given k, one may perform the singular value decomposition (SVD) of the matrix A. 
The smallest singular value is the one which we care about. If it is a minimum at a 
given k, then the billiard system is likely to have an eigenvalue at that energy. 

In the third step, the left and right eigenvectors of the smallest singular value 
($s and Cj, resp.) are used to construct a wavefunction 'J',, through a linear 
transformation. We select a grid of points Xr on which the wavefunction '^r = '^{Xr) 
is calculated. In the DEM, the left eigenvector C is used for the purpose of this 
construction, and the linear transformation which is applied is: 

*r = 5]C,F^, (6) 

3 

where Fjr = Fj{Xr \ k). Note that C contains the expansion coefficients of ^(x) in the 
chosen basis F.j{x] k). For the BIM, the right eigenvector <I> is used in order to build 
the wavefunction, and the linear transformation in this case is 

= ^G,,$„ (7) 

s 

where G^s is the discretized version of the Green function. Thus, the vector $s 
represents a 'charge' that is distributed along the boundary. 

In the final step, a measure S{k) is defined such that S{k) = if fc is an eigenvalue 
and S{k) > otherwise. In practice, the eigenvalues are determined by searching for 
the local minima of S{k). By construction, the wavefunction which was built in 
step three satisfies the Helmholtz equation inside the boundary. Therefore, the most 
natural choice of S{k) is the tension, the sum of the square of the wavefunction along 
the boundary. The tension is thus a measure for the roughness of the constructed 
\E'(a;) along the boundary. This definition of S{k) is traditionally used with the 
PWDM. Other possibilities for the measure include the smallest singular value, and 
the Fredholm determinant of A. These two latter choices of S{k) are the ones that are 
usually associated with the BIM. In Sec. 21 we discuss the mathematical equivalence 
of the three possible measures, and compare their respective numerical effectiveness. 

1.1. Outline 

In Sec. 121 we give a concise presentation of the BIM and the related GFM. Our 
derivation of the BIM equation contains some significant improvements over previous 
ones. Most importantly, it naturally leads to the existence the GFM. Furthermore, 
we have succeeded in avoiding the use of the complicated "regularized" method of 
images, which was the major ingredient in the derivation of Ref. |12j . 

Strategies for constructing the wavefunction are discussed in Sec. O An 
explanation of the Green function method is given, as well as a critical discussion 
of the DEM and its numerical variants. 
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Sec. 01 explores the practicality of using different choices for the quantization 
measure. In particular, it is demonstrated that a tension measure can be defined not 
only for the PWDM, but for the case of the BIM as well. An important issue emerges 
as to whether the quantization measures can be used as to determine the error in the 
bulk wavefunction. We address this issue, and also make a comparison between the 
numerical accuracies of the BIM and of the PWDM. 

SecEl explains how the GFM bridges between the BIM and the DEM. It is found 
that for any DEM, an associate GFM exists, whereas the inverse statement is not 
true. 

The shape that we have studied numerically is presented in Fig.l. We have used 
the cornerless, generic 'Pond' shape in order to avoid the range of problems that arise 
with more complicated geometries. These problems are the subject of a follow-up 
study 53]) where we suggest mixed BIM/DEM methods for finding eigenfunctions. 
This is done using the above theoretical framework, while regarding the 'Pond' shape 
as a reference against which to judge the effectiveness of our efforts. Another direction 
of research is related to billiards in magnetic fields |14| . 

For the convenience of the reader, our numerical notations are concentrated in 
Table 1. Further information about Fig. 1, Table 1, and the numerical analysis is 
integrated within the main text. 

2. The BIM and the GFM 

The gist of the BIM is that, from the knowledge of the gradient of the wavefunction 
on the boundary and from Green's theorem, it is possible to find the value of the 
wavefunction everywhere inside the billiard. We give a derivation of this method in 
this section. This procedure will lead us naturally to the existence of the GFM. 

The free space Green function G{x,x') is defined by the equation TCG{x,x') = 
6{x — x'). The most general solution can be written as 



where C{x,x') is any solution of the homogeneous equation Ti.C{x,x') — 0. Note that 
in the electrostatic limit fc ^ we have G{x,x') — — (l/(27r)) ln(r) +C, where C is a 
constant |2S1. We shall refer to the choice of C{x,x') as gauge freedom. This term is 
at the core of the GFM. 

By the definition of the Green function, it follows that a solution of the generalized 
Poisson-Laplace (GPL) equation H'i'ix) — p{x) is [201 



We refer to p{x) as the 'charge density', by analogy to its electrostatic equivalent. 

We shall use the notation <I>(s) in order to refer to the (surface) charge density 
upon the boundary. In the latter case, the equation above reduces to 



G{x,x') 



-Yo{k\x-x'\) + C{x,x') 



(8) 




(9) 




(10) 



where s parameterizes the boundary. 
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Xs 


— vector of bouiids-ry points 


Xs 


— vector of outer- boundary points 




— vector of interior grid points 




= randomly selected interior point 


\T/ 


— wavefunction on tlie grid points 




— wavefunction on the boundary points 
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ll*.|| 


— tension of the wavefunction (see Sec. Ill) 


n{s) 


= unit normal at the boundary point Xs 




= (l/(2fc2)) • 


Grs 


= GiXr,Xs) 






Ajs 


= Fj{xs;k) 


Dj. 


= dFj(xs]k) 


By 


= As X;.w,Di,Dj, = As (DwDt)y 




= As X^fiAisAjs = As (AAt)y 



Table 1. Notations 



^J. The BIM 

Let us assume that k is an eigenvalue of the biUiard. In such a case, there exists a 
non- vanishing ^'(x) inside the boundary that satisfies "^{x) =0 on the boundary. It 
can be shown from Green's theorem that the interior wavefunction satisfies Ea. l|lU|l 
with 

$(s) = a_*(x(s)) = hm n(s)-V*(x) (11) 

x\x{s) 

where n(s) is the outward pointing normal at point s, and d- is used for the normal 
derivative evaluated inside of the billiard walls. 

In electrostatics, it is known that forcing the scalar potential to be zero on the 
boundary induces a boundary charge. From Green's theorem, the induced charge is 
proportional to the normal component of the electric field. Here the wavefunction 
acts as the equivalent of the scalar potential. Similarly to the electrostatic case, there 
exists an induced 'boundary charge', that is in this case proportional to the normal 
derivative of the wavefunction. 

The BIM is based on the fact that if an eigenstate exists, then there also exists a 
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charge density $(s) given by Ea. l{TT|l . such that Ea.({Tn|l is satisfied. On the boundary, 
Eq.lini) yields 

J G{x{j),x{s))<i>{s)ds = [BIM equation] (12) 

Thus, having an eigenstate ^>{x) implies that the kernel G{x{j),x{s)) has an 
eigenvector $(s) that corresponds to a zero eigenvalue. Fig. 2 shows an example of a 
boundary charge density $(s) that was found via the BIM equation (for more details, 
see next section). The converse is also true: Once a non-trivial charge density is found 
that satisfies Ea. l|12|l . the associated eigenstate can be constructed using Ea. p()|l . We 
discuss this construction issue in more details in the next section. 

For numerical purposes, it is convenient to use the discretized version Eq.(|7|) of 
the above formula. The BIM equation can then be written as the matrix equation 
A$ = 0, where Ajg = G{x{j),x{s)). The gauge term C{x,x') allows some freedom 
in the determination of the Green function. Using the Neumann Bessel function 
lo(^|a^ — x'\) for the Green function, one obtains the matrix A^s as defined by Eq.(|SJl 
with lO. Another possibility is to use the Hankel Bessel function HQ{k\x — x'\). 
Accordingly, we will distinguish between the YO-BIM version and the HO-BIM version. 
We shall later discuss the numerical implication of using the complex H{k\x — x'\) 
rather than the real Y{k\x — x'\). 

The primitive BIM uses Ea. p2ll literally. However, this version of the BIM is not 
the one that is generally favored because G{x{j), x{s)) is singular for x{j) — > x(s), 
leading to some difficulty in determining the diagonal matrix elements of Aj^. 
Therefore, other versions of the BIM have become popular (sec Appendices B,C). 

2.2. The GFM 

The GFM is a different strategy to obtain the charge density $(s). Rather than using 
the BIM equation Eq. (|12|l or one of its variants, a gauge freedom argument is invoked 
in order to introduce a new type of equation fEa. H13|l below). It is clear that Ea. H12|l 
should be valid for any choice of gauge. In other words, Ea. (|12|l should be satisfied 
for any Green function (Eq.®), whatever the choice of C{x,x'). Thus, for a given 
C{x,x'), the charge density ^{s) must satisfy the equation 

J C{x{j),x{s))<^>{s)ds = [GFM equation] (13) 

For example, we may take C{x,x') = Jo(fc|a; — x'\), and we shall refer to this version 
of GFM as JO-GFM. For numerical purposes, it is once again convenient to discretize 
the integral expression, which can then be written as the matrix equation A<f> — 0. 

The kernel Ajg = Jo{k\x — x'\) of the JO-GFM is non-singular, and very well- 
behaved. Thus, the JO-GFM method, unlike the YO-BIM, provides an extremely 
convenient way of obtaining the eigenvalues of the Helmholtz equation. Fig. 2 shows an 
example of a charge density that was found via the JO-GFM equation (for more details 
see next section). The result is indistinguishable from the charge density generated 
by the traditional Hl-BIM. [We note however that the JO-GFM method has certain 
numerical limitations that we are going to discuss later]. Once the eigenvector $(s) 
is found via the GFM equation, we can proceed as with the traditional BIM, and 
construct the wavefunction "^(x) using Eq.|7I). 
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3. Constructing the wavefunction 

In this section we explain how a wavefunction '^(x) is constructed for a given k. It 
is assumed that k is an eigenvalue. The (numerical) question how to determine the 
eigenvalues k — kn is differed to Section 4. 

3.1. Green function method (Eg. ^) 

Both the BIM and GFM make use of Eq. Q in order to construct the wavefunction. In 
order to find the charge vector <I>s the BIM equation fEa. ((T^ ) and the GFM equation 
('Ea. H13ll ') are written as the matrix equation A$ = 0. The only difference between 
the two is in the expression for A. Given k, one performs the SVD of the matrix A. 
Fig. 4 displays an example of the output of such a SVD procedure. One then finds the 
right eigenvector $ that corresponds to the smallest singular value. 

Once the charge vector $s has been determined, as in the example of Fig. 2, one 
can construct the wavefunction using Eq.Q. For the Green function Eq.®, it is most 
natural to use the simplest gauge (C = 0). If fc is known to be an eigenvalue, then 
any gauge should give the same result, and in particular, the wavefunction associated 
with any complex part of the Green function (such as that of the Hankel function) 
should vanish. The outcome of the Green function method is illustrated in Fig. 3. 

It is natural to ask how the constructed wavefunction ^'(a;) look like outside of the 
boundary. The answer turns out to be ^'(x) — 0. For completeness, we give a proof 
of this statement. Let us define an extended function '^/^^{x) such that '^^^{x) = '^{x) 
inside and ^„^{x) = outside of the boundary. We would like to show that '^{x) as 
defined by Ea. ((Tr)|l is also equal to '^„^{x) outside of the boundary. It is clear that '^{x) 
is a solution of the GPL equation by construction [see discussion following Eq. . 
In the next paragraph, we are going to argue that '^„^{x) is a solution of the same 
GPL equation. It follows that the difference R{x) = ^'(x) — ^ex(a;) is a solution of 
Helmholtz equation in free space. From the definition of 'i'^^{x), we have R{x) = in 
the interior region, which implies by the unique continuation property that R{x) = 
over all space. 

The proof that 'i>^^{x) is a solution of the GPL equation with a charge density 
given by Ea. l(TT)) goes as follows: By construction, ^^„{x) satisfies the GPL equation 
inside as well as outside of the boundary. All we have to show is that it also 
satisfies the GPL equation across the boundary. The latter statement is most easily 
established by invoking Gauss' law. This approach is valid because at short distances, 
G{x,x') coincides with the electrostatic Green function. Thus, the gradient of ^'ex(a;) 
corresponds, up to a sign, to the electric field. Gauss' law implies that the electric 
field should have a discontinuity equal to the charge density $(s). Indeed, '^„^{x) is 
consistent with this requirement. 

3.2. Decomposition method (Eq. 

The other procedure to construct the wavefunction is to use the DEM Eq.®. The 
idea is to regard Fj{x; k) as a basis for the expansion: 



j 

Any such superposition at fixed fc is a solution of the Helmholtz equation within the 
interior region. Thus, in order to satisfy the Dirichlct boundary conditions, one looks 




(14) 
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for a vector C of expansion coefficients that satisfy CA = 0. It turns out that the direct 
numerical implementation of this simple idea is a complicated issue (see discussion of 
the null space problem later in this section). 

Any set of Fj(x]k) which are solutions of the Helmholtz equation may be used 
for the DEM. However, it should be remembered that computationally not all bases 
are equivalent. For instance, the Yq basis defined by Eq.Q), which might appear 
to be the best choice as a DEM basis due to its association with the BIM, does not 
give the best numerical results when compared against other options. In particular, 
it turns out that the PWDM is typically much more effective (recall that the PWDM 
is a special case of the DEM, corresponding to the choice Eq.(|2l of basis functions.) 
Finally we note that the set of JqS of Eq.Q) can not be regarded as a mathematically 
legitimate basis for a DEM. This latter point will be explained in Section 5. 

For the Yq basis the BIM and the DEM lead to the same equation. It is only the 
mathematical interpretations that is different. Within the DEM, one regards the Yq 
as basis functions to be used in an expansion, while the same Yq in the BIM context 
serves as the Green function. In the context of DEM, one may be bothered by the 
singular nature of the Yq functions: The constructed wavefunction should be zero 
on the boundary. Mathematically this is achieved in the N ^ oo limit, so the Yq 
basis is a valid choice. But in an actual numerical implementation, the wavefunction 
so constructed will always have singularities on the boundary. One possible remedy 
consists of enforcing the boundary conditions on intermediate boundary points, or on 
points that lie on an outer boundary. Alternatively, one may replace the bare Yq basis 
by smooth superpositions of Iq functions (see Appendix C). 

In Fig. 2, we display an example of a numerically determined C (using the PWDM 
basis) for one of the Pond eigenstates, while in Fig. 3 we illustrate the constructed 
wavefunction. Unlike the Green function construction, the DEM/PWDM constructed 
wavefunction does not vanish outside of the boundary. Actually, it is quite the 
opposite: Typically the DEM/PWDM wavefunction becomes exponentially large as 
we go further away from the boundary. Whenever this behaviors occurs, it constitutes 
an indication of the evanescent nature of the wavefunction. Namely, in such cases, the 
wavefunction acquires sub- wavelength features in order to accommodate the boundary. 
This requires exponential behavior (negative kinetic energy) in the transverse space 
direction, in order to keep the total energy fixed. 

3. 3. Normalization of the wavefunction 

A standard SVD procedure generates vectors <I>s and Cj that arc normalized in the 
sense \^s\'^ — 1 and |Cjp — 1. Therefore, the constructed ^'r is not properly 
normalized within the interior region. Adopting the usual philosophy of boundary 
methods, the problem of calculating the \E'r normalization is reduced to that of 
evaluating a boundary integral, namely |H1 



For the BIM, by discretizing of Eq. H15|) we obtain the following numerical expression 
for the normalization factor: 




(15) 




(16) 



S 
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Here W — (1/As)diag(ws) is a diagonal matrix, and the weight factor Wg is defined in 
Table 1. As for the DEM, by using the derivative of Ea. ((TH) in Ea. lfTTTl and substituting 
into Ea. ((T^ . we get: 

= E^^B.^C. = CBCt (18) 

The definitions of Djs and of the metric Bij can be found in Table 1. 

The normalization jlvPrll can be calculated using the metric By. This method 
looks quite elegant, but it turns out not to be very effective numerically. Consider 
Eg. ljlTI) . This equation is quite safe computationally for two reasons: (i) all its terms 
are non-negative; (ii) standard summation routines order these terms in descending 
order. Now, let us look instead at Ea. H18|l . In this case the numerical calculation can 
give any result (if we go to large k). Sometimes, the answer even comes out to be 
negative! This occurs because the calculation involves many arbitrarily ordered terms 
that each have a different algebraic sign. 



3.4- The tension along the boundary 

The numerical wavefunction \E'r, satisfies Helmholtz equation in the interior region by 
construction. Thus, whether is an actual eigenstate depends on its behavior along 
the boundary. In this subsection we would like to discuss the definition of a 'tension' 
measure that estimates whether, and to what accuracy, the numerical satisfies the 
boundary conditions. 

For the case of the DEM, following j9j, the tension is defined as the boundary 
integral 

||vl/,|| = AsJ] j^QA.J = 5]QT,,g = CTCt (19) 

s \ j / ij 

The standard practice to date for the tension calculation has been to use a denser 
set of boundary points located between the a;^ points. However, our experience (see 
also 0) is that the tension estimate obtained from the initial set of points is just as 
effective. This is demonstrated in Fig. 3d. Therefore, we routinely rely on the same 
set of boundary points to determine the tension. 

For the BIM on the other hand, the above definition is not practical due to the 
singular nature of the basis functions. For any finite M , the numerical wavefunction 
blows up at each boundary point. However, since the BIM wavefunction vanishes 
everywhere outside of the billiard, a numerically unambiguous definition of tension 
arises as an integral of |^(a;)|^ along an outer boundary: 

W^sW = As ^(vI/,)2 = A,s ^ (E^^"*^) (20) 

By outer boundary (see Fig.l), we mean the set of external points (s points, as opposed 
to s points for the true boundary) that have a fixed transverse distance Ai from the the 
true boundary. The distance AL between the boundary and the outer boundary should 
be small on any classical scale but large compared with ds, in order for the tension to 
be independent of the choice of AL. See Fig. 5c for a numerical demonstration. 
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3.5. The PWDM and the null space problem 

One may think that C could be found simply by computing the (left) eigenvector 
that corresponds to the smallest singular value of Ajg. Numerically this definition is 
hard to implement. This difficulty can be explained by looking at the behavior of the 
singular values of Ajs for the PWDM basis. Fig. 4 gives some examples of singular 
values resulting from the SVD of the Ajs matrix. In the case of the PWDM, as 
k become large, one observes that the the singular values separate into two groups: 
rather than having one distinctly smaller singular value, we obtain a whole set of them. 
Accordingly, we can define a numerical 'null space' of the Ajs matrix. 

The interpretation of this null space is quite clear. It is well known |19[ I20j that 
it is not efficient to include much more than Ngc plane waves in the basis set Fj{x; k), 
where 

Nsc = -kL (21) 

TT 

and L is the perimeter of the billiard. The reason for this ineffectiveness is that ki and 
kj cannot be distinguished numerically within the interior region unless \ki — kj\L > 1. 
In order to obtain the semiclassical result (|21|) . the total phase space area (L x {2mv)) 
of the boundary Poincare section is divided by the size of Planck cell {2'itK). If we use 
N > Nsc plane waves, then we can create wavefunctions that are nearly zero in the 
interior, and become large only as we go far away from the center (21). It is clear that 
the SVD can be used to determine the N — Nsc null space of these evanescent states. 
Whenever k is an eigenvalue, this null space includes, in addition to the evanescent 
waves, the single eigenvector that constitutes an eigenstate of the Helmholtz equation. 
The problem is to distinguish this eigenvector from the other vectors in the null-space. 

The basic difference between the eigenvector that leads to an eigenstate (which 
exists if fc = kn), and the other vectors in the null space is related to the normalization. 
As discussed before, a standard SVD procedure generates vectors Qj that are 
normalized in the sense |Cj|^ = 1. Therefore, the ^'^ of Eq.ljSJ is not properly 
normalized within the interior region. Normalizing the wavefunction has the effect of 
magnifying the evanescent solutions in the interior as well as on the boundary, while 
the eigenfunction (if it exists) remains small on the boundary. In appendix D, we give 
a detailed explanation of the numerical procedure for finding Cj that can be derived 
from the above observation. 



4. The quantization measure 

Once we have constructed the wavefunction at a given fc, the next step is to determine 
whether 5* is an eigenstate. As we will explain below, our choice of measures reduces 
to finding the minima of one of: 

S{k) = tension (22) 
S{k) = smallest singular value (23) 
S{k) = determinant (24) 

We call S{k) the quantization measure. Below we give further explanation of the 
above definitions. 

It is clear that the most natural quantization measure is the tension. If a properly 
normalized wavefunction has "zero tension" on the boundary, it means that the 
corresponding k is an eigenvalue. The normalization issue is further discussed in 
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Appendix D. The question that arises is whether we can use a numerically simpler 
measure, and what price we pay for doing so. 

The BIM Eq. ^ and the GFM Eq. (O can both be written as A$ = 0, with 
the appropriate choice of A. Thus, if k is an eigenvalue, A should have a singular value 
that tends to as A'' increases. The determinant of A is defined as the product of all 
the singular values, and therefore it should vanish whenever one of the singular values 
does. Using the GFM-DEM duality which is discussed in Section 5, it is clear that for 
the DEM (and for the PWDM in particular) the determinant of A vanishes whenever 
k is an eigenvalue. It is important to realize that in the latter argumentation, we do 
not rely on inside- outside duality |20| considerations, but rather on the much simpler 
GFM-DEM duality. 

Thus, a low tension must be correlated with having a vanishingly small singular 
value or determinant. The converse is not true: It is well known that SVD based 
quantization measures may lead to spurious minima (see 0] and references therein). 
Therefore SVD based procedure for finding eigenvalues requires a post-selection 
procedure whose aim is to distinguish true zeros from fake ones. 

It is important to realize that neither the traditional implementation of PWDM, 
nor that of the BIM should be considered to be 'package deals'. For example, the BIM 
could be used with the tension as a measure (defined in the next section) , rather than 
looking for minima of the the singular values. Similarly, the usual Heller method 
of PWDM implementation (see Appendix D) could be replaced by a search over 
determinant values. 

4-1. The tension as a quantization measure 

The tension is a robust measure of quantization. Fig. 5 displays some examples of the 
corresponding S{k) plots. The PWDM minima are typically much sharper than their 
BIM equivalents. Zooming over a PWDM minimum (Fig.5d) reveals some amount 
of roughness. This feature is actually helpful, because it gives an indication of and 
control over the accuracy of the numerics. We interpret the roughness of the PWDM 
minimum as a reflection for the existence of a null-space. In the same spirit, the 
smoothness of the BIM minima can be regarded as an indication that better accuracy 
can be obtained by making N larger. We discuss these issues below. 

The tension provides a common measure that may be used to monitor 
improvements, as well as to compare the success of the different methods. Naturally, 
the first issue to discuss is the dependence of the tension on the size N of the basis set 
(see Fig. 6). For the BIM, the tension becomes better as N grows, and disregarding the 
computer hardware, there is no reason to suspect that there is an inherent limitation 
on the accuracy. The situation is different for the PWDM. Here, taking TV much larger 
than Nsc is not effective. In practice, the method reaches a limiting accuracy, which, 
taking into account present hardware limitations, is still very good compared with 
that of the BIM. 

From Fig. 6, it is also clear that the tension of the PWDM becomes much better 
as k becomes larger. This is expected on the basis of the following semiclassical 
reasoning: larger uncertainties in k result from confining a particle to a smaller box 
(taking a smaller box for a given k is equivalent to making k smaller for a given box 
size). Thus, it is more difhcult to build a wavefunction with a precise value of \kj \ = k 
for low lying eigenstates. On the other hand, the BIM does not seem to be sensitive 
to the value of k. 
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4-2. The tension as an indication for the global error 

The tension can be regarded as a measure of the local error in the determination of 
the cigcnfunction. The tension is local in the sense that it pertains only to points 
along boundary. We can also define a measure for the global error, that is the error 
which is associated with all the interior points: 

(A*)2 = (25) 

Here ^'exct(a^) is the numerically exact wavefunction. The average is taken over the set 
Xr of selected points inside of the boundary. Fig. 7 gives an example for the variation 
of the error along the cross section line of Fig.l. In order to eliminate a possible bias 
due to a global normalization error, we renormalized the inexact wavefunction so that 
= '^c,,^t{Xr) at a randomly selected point X = Xq. In retrospect, we realized that 
such an error did not significantly affect the result. However, we still chose to be on 
the safe side, and we adopted this procedure routinely. 

It is natural to expect the average error (A^*)^ to be correlated with the tension. 
In other words, if |\l/r- — ^'exct(-^r)| is small on the boundary, then one may expect it 
to be small in the interior. The degree of such correlation is important for practical 
reasons. Moreover, we have introduced two different versions of tension definitions, 
one for each of the PWDM and the BIM. It is not a-priori clear that the above 
correlation is independent of the choice of numerical method. In Fig. 8, we study this 
issue by plotting (A^*)^ against the tension for the BIM and the PWDM. In the case 
of the PWDM, the error saturates below a critical tension. After this point, further 
improvements on the boundary do not seem to affect the bulk of the eigenstate. It 
is not clear from the numerics whether or not the BIM saturates. What is clear 
however is that the BIM does a poorer job at reproducing the wavefunction inside of 
the boundary than the PWDM with the same tension. 

The satiiration of the error well inside the billiard can be explained as a 
manifestation of the fact that the wavefunction there is not very sensitive to sub- 
wavelength roughness of the boundary: If N is reasonably large, the numerical 
wavefunction vanishes on a nodal line that almost coincides with the true (pre-defined) 
boundary. Increasing N further affects the sub-wavelength features of the (distance) 
difference between that nodal-line and the true boundary. This distance difference is 
important for the tension, but barely affects the wavefunction well inside the billiard. 

4-3. The determinant as a quantization measure 

The tension is the natural choice for a quantization measure. However, from a 
numerical point of view, it is much more convenient and time effective to compute 
the singular values of A, without having to find the eigenvectors for each k value, 
and without having to compute the wavefunction along the boundary (for tension 
calculation) . 

The smallest singular value is traditionally used as a quantization measure for 
the BIM. Prom Fig. 7, it is quite clear that for the BIM one of the singular values 
is significantly smaller than the others, so that the eigenstate is unambiguously 
determined by this method. Is it possible to use the same approach, with comparable 
success, for the PWDM? We have already determined that looking at the smallest 
singular values is not very meaningful numerically. For A'^ > TV^c, there exists a 
large null-space of evanescent states for any k. The metric method (Appendix D) is 
not a practical solution to this problem since we want to gain numerical efficiency 
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(if efficiency is not the issue then it is better to use the tension as a quantization 
measure). 

One simple way to improve the numerical stability is to use the determinant rather 
than the smallest singular value as a quantization measure: Each time that k = kn, the 
null space should include one more 'dimension'. Therefore, the determinant, rather 
than the smallest singular value, becomes the reasonable quantity to look at. Thus, 
from numerical point of view 1)23(1 should be superior compared with ((22|l . 

Fig. 9 illustrates how the determinant can be used in practice as a quantization 
measure. As a general rule, as is the case for the tension, the PWDM/GFM minima 
are sharper than the BIM ones. On the one hand, this extra sharpness can be regarded 
as an advantage, because it leads to a better resolution of the eigenvalue spectrum. 
However, more computer time is needed in order to find these minima. The BIM 
minima are broader, and therefore digging algorithms that search for local minima 
become extremely effective. 

In the case of the traditional BIM, using a larger N leads to a better resolution of 
the local minima, as expected. The traditional H-BIM uses the complex Hankel Bcssel 
fimction as its Green function. One may wonder why the real Neumann function could 
not be used instead. A-priori, there is no reason to insist on Hankel choice. However, 
it seems that with Neumann choice the numerics are not very stable: The locations of 
the local minima vary on a fc range which is large compared to their k width. Because 
of this problem, search routines relying on the Y-BIM may yield misleading values for 
the error in the fc„ determination. Thus, the numerical stability of the H-BIM can be 
attributed to the fact that the BIM equation A$ = becomes complex. Its real part 
is just the Y-BIM equation, while its imaginary part is the J-GFM equation. Thus 
one may say that the H-BIM benefits from combining the Y-BIM with the J-GFM. 

Is it practical to use the Fredholm determinant as a quantization measure also in 
the GFM/PWDM case? Here we observe that the null-space problem is reflected in 
the stability of the determinant calculation. It is useful to characterize the numerics 
using the discretization parameter b: 



The last equality holds if we take M — N, leading to the interpretation of h as the 
number of boundary points per half De-Broglie wavelength. If 6 < 1 the null space 
problem does not exist, and we can define C as the (left) eigenvector that corresponds 
to the smallest singular value of Ajg. Of course, we want to push PWDM to the limit, 
and therefore in practice we always take b > 1. 

The natural question is whether choosing a very large b is numerically useful. 
Our numerical experience is that for 1 < 6 < 1.8 we get nice minima, which actually 
look much sharper than the BIM ones (see Fig. 9). As we try to increase b in order 
to improve accuracy, the numerics loose stability (what we mean by instability is 
demonstrated in Fig. 9e). The same phenomenon occurs with JO-GFM, which has 
somewhat larger tendency for instability. This is apparently because the JO-GFM is 
involved with a larger null-space (see Fig. 4). 

Thus we conclude that taking 6 > 1 does improve the numerics, while 6 3> 1 
generally leads to instabilities that should be avoided. The optimal choice of b 
depends on the details of the implementation and on the computer hardware. It 
should be clear that if the numerical accuracy were unlimited, then the 6 — > oo 
limit would lead to a numerically exact solution in cases where the wavefunctions 




A/2 



(26) 



M = N 



As 
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may be written as superpositions of plane waves. This is not always possible jl7j . 
Note however that evanescent features of the wavefunction can be reconstructed by a 
suitable superposition of plane waves |21| . 

5. The duality of the GFM and the DEM 

In addition to providing a boundary method of its own, the GFM also serves to bridge 
the gap between the BIM and the DEM. Consider the version of the GFM that is based 
on the choice C{x, x') — Fj{x; k), where the Fj are solutions of the Helmholtz equation 
in free space (with neither singularities nor cuts). With this choice, we immediately 
realize that the DEM and the GFM are dual methods: 

A$ = [GFM equation] (27) 

CA = [DEM equation] (28) 

The only difference lies in whether one looks for the left or the right eigenvector. This 
point is numerically demonstrated in Fig. 2. 

The PWDM version of the DEM also satisfies this duality. In this special case, 
a somewhat more elegant version of the above argument is as follows: Consider the 
version of the GFM that is based on the choice C{x, x') — exp{iknj ■ (x — x')), where 
rij is a unit vector in a given direction. We can take N different choices of Uj, thus 
obtaining the matrix equation = with 

Ajs — exp^ikuj ■ Xs) (29) 

An equivalent matrix equation is found by multiplying each equation by exp(i/c0j), 
where are random phases. We can then take the real part of these equations, 
thus obtaining a set of equations that involves the same matrix A as that of PWDM, 
namely © with the basis defined by Q. 

The duality of the PWDM and the GFM is very important from a mathematical 
point of view. The mathematical foundations of the PWDM are quite shaky. It is 
clear that PWDM is well-established mathematically if a strict 'inside-outside duality' 
(lOD) 1201 is satisfied. In this case, the wavefunction can be extended to the whole 
plane so that the boundary of the billiard can be regarded as a nodal line of some 
plane- wave superposition. Obviously, this is rarely possible [21] ■ Therefore, one may 
wonder whether we indeed have det(A) — whenever k = kn- Using the duality 
Ea. H27|) . it becomes obvious that the Fredholm determinant indeed vanishes at the 
eigenstates, even in the absence of an exact lOD. 

It is quite clear from the first paragraph of this section that any expansion method 
can be associated with a corresponding GFM. Whenever the left eigenvector is used 
with the expansion method, the right eigenvector can be used with the GFM. We have 
already demonstrated this point in Fig. 2. Is it possible to make the inverse statement? 
Do we have a well defined expansion method associated with any GFM? The answer 
is negative. We discuss this issue in the rest of this section, and it can be skipped at 
first reading. 

For the following, it is convenient to consider Eq.® as — > oo. Subsequently, 
we are going to talk about whether this limit is meaningful. In the case of the usual 
PWDM, the N oo limit of Eq. I© can be written as 

*(a;) = / C{e)de exp{ikn{e) ■ x) (30) 
Jo 
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Similarly, in the case of the Jq decomposition, using the basis functions of Eq.lQ, we 
can write in complete analogy: 

^(x) = (j) ^{s)ds Jo{k\x - x{s)\) (31) 

In writing Ea. H31|l we have used the fact that Jo{x{j) — x{s)) is a symmetric kernel, 
and therefore we could make the substitution C = 

Eq.JSU looks at first sight like an innocent variation of Ea. H30|) . The Bessel 
function Jo is just a superposition of plane waves, and therefore one may possess the 
(incorrect) idea that there is a simple way to go from Ea. (|31|l to Ea. (|30|l . If we expand 
each Jo in Ea. pil) in plane waves, and re-arrange the expression in order to identify 
the PWDM coefhcients, we end up with the relation 

C{0) = J e-*'^'"(»)-^(^) $(s)ds (32) 

This relation implies the trivial result C{9) =0 due to gauge freedom [see discussion 
of Ea. (j29|l ]. Hence we conclude that the constructed wavefunction is 4'(x) =0 in the 
N —f oo limit! 

Having "^(x) = from Ea. H31() could have been anticipated using a simpler 
argument: We know that Eq.l^UJ should hold for any gauge choice. This gauge 
freedom implies that we have 

Cix,x{s))<P{s)ds = (33) 

The above equation should hold for any x inside as well as on the boundary. 
Furthermore, the left hand side of 1)33(1 is manifestly a solution of Helmholtz equation 
in free space, and it follows from the unique continuation property that Eq. I(33() holds 
also for points outside of the boundary. Having ^'(a;) = as a result of the integration 
in Ea. (|^ is just a particular case of Ea. llH^ . 

In spite of the observation that Ea. (|31|l yields 'i>{x) = in the TV — > oo limit, 
the vector is non-zero numerically for any finite N. In fact, after proper re- 
normalization, "ifr becomes a quite good approximation to the wavefunction (see Fig. 3c 
and Fig.3e). Sometimes the result so obtained is even better than the one which is 
found via the traditional BIM Eq. |7J). As strange as it sounds, this success is entirely 
due to the fact that we are using finite iV. 



6. Conclusion 



The BIM and the DEM were written as different faces of a unified boundary procedure 
comprising four steps. In the process of doing so, yet a third boundary method was 
derived as part of the same framework, the GFM. The DEM and the GFM are strongly 
related, as they are respectively the left and right eigenvectors of the same Fredholm 
matrix. We think that the presented approach opens the way towards a controlled 
fusion of the BIM and the DEM into a more powerful numerical procedure. 

The unified treatment of quantization measures allowed us to compare the 
efficiency of the various methods, and to analyze both the local and the global errors in 
the numerically determined wavefunction. In particular, a numerically valid definition 
of the BIM tension was given, and was found to possess smooth minima at the 
eigenstates. Using the tension as a quantization measure is one possible way to avoid 
some problems 0] that are encountered in the traditional implementation of the BIM. 
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Appendix A. The BIM for the scattering problem 

The solution of the Helmhohz equation for the scattering problem is just another 
variation of the BIM. Consider a boundary, one that in general may be composed of 
several disconnected pieces. The incident wave *I'i„cidont(a;) is a solution of Hclmholtz 
equation in free space. Formally, ^^incidont (2:) includes both the ingoing and the outgoing 
wave components. We look for a solution \E'(x) that has the same ingoing component 
as 4'i„cidont(a;): and that satisfies ^'(x) = on the boundary. Such solution can be 
written as a sum of the incident wave and a scattered wave, and hence must be of the 
form 



Equation ljA.l|l is a variation of Eq. . Note that the Green function should satisfy 
outgoing boundary conditions in order to yield the desired solution. The charge density 
$(s) is fixed by the requirement that ^'(x) = 0, which leads to the boundary equation 



This inhomogcncous equation is a straightforward generalization of Ea. H12|l . A 
discretized version of it was used in Ref . j22j in order to obtain numerical solutions 
of the Helmholtz equation for some scattering problems. 

The derivation of Ea. ljA.2|l in Ref.[22| is much more complicated than ours, and 
involves the use of Lippmann-Schwinger equation. The boundary is represented by 
a large delta-potential V, and the limit ^ ^ 00 is taken. Using this procedure, the 
charge density <i>(s) can be obtained as the V ^ 00 limit of V'i/{x{s)). Note the 
correctness of the physical units. Namely, = [p] and therefore = [<&]. 

Note also that the wavefunction is in general non-zero in both sides of the boundary. 
Therefore the charge density <&(s) is equal to the difference between the normal 
derivatives on both sides of the boundary. The simplest way to derive the relation 
between $(s) and V'^{x{s)) is to integrate the Helmholtz equation over an infinitesimal 
range across the boundary, as in the treatment of the ID Shchrodinger equation with 
delta potential. 

Appendix B. Traditional BIM equation 

The primitive BIM equation (Eq.l^SJ) is based on Eq.l^Hl- The traditional BIM is 
a variation of the same idea which avoids the boundary singularities that plague the 
more simplistic version. Rather than using Eg. ljlUI) directly, one considers its gradient, 
leading to 



This equation is analogous to Ea. (|12|l . We use the notation 9+ and 9_ in order to 
refer to the normal derivative 

on the interior and exterior sides of the boundary By definition, 9_5'(x(s)) = ^{s), 
and from the discussion in section 111 we have (?+5'(x(s)) = 0. Adding the two equa- 




(A.1) 




(A.2) 




(B.l) 



tions of IjB.ip . we obtain 




(B.2) 
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where d = + 9_)/2 is just the derivative on the boundary in the principal sense. 
Thus, in the traditional BIM, the definition of the Fredholm matrix is 

A,,. = -^Sss' - 2dG{x{s),xis')) (B.3) 
As 



and the BIM equation (|B.2p is A$ = 0. The kernel dG{x{s), x{s')) is well- behaved, 
and its diagonal elements are finite thanks to the presence of a geometrical factor. 
Namely, if G{x{s),x{s')) = g{k\x{s) — x{s')\) then 

dG^k ' ' , ' " g'{k x{s)-xis') ) B.4 

\x{s)-x{s')\ 

If either one of the Bessel functions Hq or Yq is used for the Green function, then the 
definition of A^^/ above involves either Hi or Yi, respectively. 



Appendix C. Transformed BIM equations 

There exists another elegant version of the BIM 1151 III)] which does not exhibit 
the singularities associated with the YO-BIM. Namely, the BIM equation is rewritten 
as 

J G{x{j),K)^{K)dK = 0, (C.l) 

where the vector is related by a linear transformation to the vector ^(s). This 
transformation corresponds to the selection of a complete basis set of boundary 
wavefunctions. The KKR method |16 | is a particular implementation which uses 
the spherical harmonics. Another (more general) choice E| consists of taking 
as the Fourier components of ^{s). In the latter case G is related to G by the 
Fourier transform s i— > k. However, obtaining G from G is not a simple matter for 
most boundary shapes. 



Appendix D. The PWDM, in practice 

The mathematically clean solution for the null-space problem is to adopt a method 
based on a metric. As we explain in the next paragraph, this procedure is sensitive to 
cumulative numerical errors. A modified implementation of the metric method, that 
avoids some of the numerical problems, has been introduced by Barnett ^5, . The other 
possibility is to use a very simple procedure which is known as Heller's method jS]. 
Below we discuss the latter as well. 

The metric method works as follows: First, one finds the basis in which the 
normalization metric B^- becomes 6ij . The tension metric Tij should then be written 
in that same basis. The SVD is done on the transformed tension metric. In this case, 
the null space becomes at most one-dimensional (whenever k = fc„). Unfortunately, 
this elegant and straightforward metric scheme does not work very well, due to the 
finite precision problems discussed in connection with Ea. H18|l . Furthermore, B and 
T are "squares" of A, which leads to a loss of numerical precision compared with an 
A-based strategy. 

The most widely used A-based strategy is referred to as Heller's method 9 . The 
idea is to find Cj as the solution of the M > N set of equations J^j ^j-^js — 0, with 
the additional constraint J2j Q-^iO = 1- Table 1 gives the definition of Ajq. 
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By constraining the wavefunction to be '^{Xq) = 1 at a selected point Xq in the 
interior of the bilHard, we ehminate the problems associated with evanescent states, 
as the associated (evanescent) wavefunctions no longer vanish on the boundary and 
thus, there is no longer a null-space problem. As a result, quite large b can be used 
without encountering numerical instabilities. The only worry with this method is that 
Xo may happen to be very close to a nodal line. In such cases, the tension will be 
large due to an improper normalization, so we will miss these eigenstates. 
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FIG.l: Full line: The 'Pond' shape '24:' ■ Its radius is roughly 1, and its perimeter is L = 6.364. 
Dashed line: The cross section line that is used e.g. in Fig. 3c. Dots: The 'outer boundary' which 
is used for the BIM tension calculation (see Sec.III-A). The transverse distance between the actual 
boundary and the outer boundary was chosen in most cases to be AL ~ lOAs. Star: The selected 
point which is used in Heller's implementation of the PWDM method. 
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FIG. 2: The eigenvectors of the Fredholm matrix (Eq. jFJ) are found for k = kn = 
6.82754592867694. Right plot: The right-eigenvector The symbols (x) and (+) and (o) 

correspond respectively to the PWDM choice of Eq.jJJ, to the Hl-BIM choice of EQ. iB.3l . and 
to the JO-GEM choice of Eq.|4j. Left plot: The left-eigenvector C for the PWDM Fredholm matrix. 




FIG.Sabc: See captions in the next page. 
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FIG. 3: The eigenfunction at fc = fc„ is found using tiie eigenvectors of Fig. 2. (a) An image of 
'^(x) using Eq.Q. (b) The same using PWDM and Eq.JSJ. (c) Plot of ^(x) along the crossection line 
of Fig.l. The vertical lines indicate the location of the boundary. The wavefunctions that correspond 
to images a-b are shown with (+) and (x) respectively. We also show with (o) the wavefunction that 
is obtained using JO-GFM and Eq.jSJ. Panels d-e-f display plots of log{|>I'{a;)p) along the boundary. 
The symbols are as in c. 
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FIG.4: Singular values of the Fredholm matrix for kn = 6.82754592867694 (left panel) and 
for kn = 50.05474912004408 (right panel). The various symbols are as in Fig. 2. 
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FIG. 5: The tension as a function of k in the vicinity of k„ = 10.14707971517264 (panels 
(a),(c),(d)), and of fc„ = 50.05474912004408 (pancl(b)). In the upper panels (a-b) the full line is for 
the PWDM constructed wavefunction, while the dashed line is for the BIM constructed wavcfunction. 
For the low k state we chose 6 = 4, while for the high k we used 6 = 2. Panel (c) demonstrates the 
dependence of the BIM tension on the choice of the distance AL. The different curves (from the upper 
to lower) correspond to AL/As = 1, 8, 16, 12, 4. Panel (d) is a zoom over the PWDM minimum. 
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FIG. 6: The tension for the constructed cigcnstate versus the number A'^ of basis functions. 
The left panel is for the kn = 2.40425657792391 eigenstate, and the right panel is for the 
kn = 6.82754592867694 eigenstate. The symbols (o) and (+) are for the BIM and for the PWDM, 
respectively. 
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FIG. 7: Plot of the error |*,. - 'iexct{Xr)\^, 
along the cross-section line of Fig.l. We 
refer here to the A,„ = 6.82754592867694 
eigenfunction. The numerically 'exact' 
wavefunction is our best PWDM-constructed 
wavefunction (TV = 69) with tension= 10~^^. 
The 'non-exact' wavefunction is either 
BIM-constructed (-I-) or PWDM-constructed 
(x), with tension ^ 10~®. In the middle 
point the error is zero by construction (see 
explanation in the text). 



FIG. 8: The averaged error in the deter- 
mination of the wavefunction, versus the 
tension for the fc„ = 6.82754592867694 
eigenfunction. For the BIM-constructed 
wavefunction we use squares and (o), while 
for the PWDM one we use (-I-) and (x). The 
error has been determined with respect to 
the 'exact' wavefunction '^exct- The latter is 
numerically defined as either the best BIM 
eigenstate (squares and (-|-)) or as the best 
PWDM one ((o) and (x)). For both choices 
*exct had a tension better than 10"^". 
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FIG.9: The Predholm determinant (S(fe)) versus k in the vicinity of fe„ = 10.14707971517264. 
Note that S{k) is normalized such that S{k) = 1 away from the minima. Panels a-b-c show S{k) in 
the cases of the Hl-BIM, the Yl-BIM and the PWDM, respectively. The lines plotted, in order of 
decreasing S{k) minimum, correspond to 6 = 2,3,4,8 in panel (a), 6 = 4, 8, 13, 12 in panel (b) and 
6 = 2 in panel (c). The PWDM run in panel (c) was repeated 3 times with different values of the 
randomly chosen plane wave phases. Panels (d) and (e) give a zoom over the minima of panels (a) 
and (c), respectively. We witness some numerical instabilities for both the Yl-BIM and the PWDM, 
though in the latter case it is much much weaker, and can be resolved only in the zoomed plot 
(panel (e)). For larger b values, the PWDM instability is enhanced, and the results are reduced to 
numerical garbage (not shown). 



